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Abstract 

Over the billions of years since the Big Bang, the lives, deaths and afterlives of 
stars have enriched the Universe in the heavy elements that make up so much of 
ourselves and our world. This review summarizes the methods used to evolve these 
nuclear abundances within astrophysical simulations. These methods fall into 2 cat- 
egories; evolution via rate equations and via equilibria. Because the rate equations 
in nucleosynthetic applications involve a wide range of timescales, implicit methods 
have proven mandatory, leading to the need to solve matrix equations. Efforts to 
improve the performance of such rate equation methods are focused on efficient so- 
lution of these matrix equations, in particular by making best use of the sparseness 
of these matrices, and finding methods that require less frequent matrix solutions. 
Recent work to produce hybrid schemes which use local equilibria to reduce the 
computational cost of the rate equations is also discussed. Such schemes offer sig- 
nificant improvements in the speed of reaction networks and are accurate under 
circumstances where calculations which assume complete equilibrium fail. 
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1 Introduction 



Research by Helmholtz, Kelvin and others throughout the second half of the 
19th century made it clear that neither gravity nor any other then known 
energy source could account for the geologically determined age of the Sun 
and solar system. Enlightened by Rutherford's 1911 discovery of the atomic 
nucleus, Eddington and others suggested that nuclear transmutations might 
be the remedy to this quandary. With the burgeoning knowledge of the prop- 
erties of nuclei and nuclear reactions in the 1930s, 1940s and 1950s came a 
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growing understanding of the role that individual nuclear reaction played in 
the synthesis of the elements. In 1957, Burbidge, Burbidge, Fowler & Hoyle [1] 
and Cameron [2] wove these threads into a cohesive theory of nucleosynthesis, 
and demonstrated how the solar isotopic abundances bore the fingerprints of 
their astrophysical origins. Today, investigations refine our answers to these 
same two questions; how are the elements we see on Earth and throughout 
the universe formed, and how do these nuclear transmutations, and the energy 
they release, affect their astrophysical hosts? 

In this article, we will review the numerical methods used to study these 
questions. Two basic numerical methods are at the heart of the study of ther- 
monuclear kinetics in astrophysics, the tracking of nuclear transmutations via 
rate equations and via equilibria. We will also briefly discuss work which seeks 
to meld these methods together in order to overcome the limitations of each. 
Additionally we will discuss the issues encountered when coupling the nuclear 
evolution to hydrodynamic simulations. 



2 Thermonuclear Reaction Networks 

Composed of a system of first order differential equations, the nuclear reac- 
tion network has sink and source terms representing each of the many nuclear 
reactions involved. Prior to discussing the numerical difficulties posed by the 
nuclear reaction network, it is necessary to understand the sets of equations we 
are attempting to solve. To this end, we present a brief overview of the ther- 
monuclear reaction rates of interest and how these rates are assembled into the 
differential equations that must ultimately be solved. For more detailed infor- 
mation, we refer the reader to several textbooks covering this subjects [3,4,5]. 
We will end this section by briefly discussing the coupling of nucleosynthesis 
with hydrodynamics. 

2. 1 Thermonuclear Reaction Rates 

There are a large number of types of nuclear reactions which are of astrophys- 
ical interest. In addition to the emission or absorption of nuclei and nucleons, 
nuclear reactions can involve the emission or absorption of photons (7-rays) 
and leptons (electrons, neutrinos, and their anti-particles). As a result, nu- 
clear reactions involve three of the four fundamental forces, the nuclear strong, 
electromagnetic and nuclear weak forces. Reactions involving leptons (termed 
weak interactions) proceed much more slowly than those involving only nucle- 
ons and photons; however, these reactions are important because only weak 
interactions can change the global ratio of protons to neutrons. 



2 



The most basic piece of information about any nuclear reaction is the nuclear 
cross section. The cross section for a reaction between target j and projectile 
k is defined by 

number of reactions target _1 sec _1 r/rij 
flux of incoming projectiles n k v 



The second equality holds when the relative velocity between targets of num- 
ber density rij and projectiles of number density is constant and has the 
value v. Then r, the number of reactions per cm 3 and sec, can be expressed 
as r = avrijUk- More generally, the targets and projectiles have distributions 
of velocities, in which case r is given by 

rj,k = / <r(\vj - v k \)\vj - v k \d 3 njd 3 n k . (2) 



The evaluation of this integral depends on the types of particles and distri- 
butions which are involved. For nuclei j and k in an astrophysical plasma, 
Maxwell-Boltzmann statistics generally apply; thus, 

rf 3 n = n( ^_)3/2 ( _^l )d 3 (3) 

y 2nk B T ! PV 2k B T J ' v ; 



allowing rij and n k to be moved outside of the integral. Eq. 2 can then be writ- 
ten as rj^k = ( av )j k n i n ki where (av) is the velocity integrated cross section. 
Equivalently, one can express the reaction rate in terms of a mean lifetime of 
particle j against destruction by particle k, 

nU) = — W 

(<Tv) jk n k 



For thermonuclear reactions, these integrated cross sections have the form [3,6] 

oo 

(j,k) = (av) jk = (A)1/2 (A . bT )-3/2 f Ea(E)eM-E/k B T)dE, (5) 



where fj, denotes the reduced mass of the target-projectile system, E the center 
of mass energy, T the temperature and k B is Boltzmann's constant. 

Experimental measurements and theoretical predictions for these reaction 
rates provide the data input necessary to study astrophysical thermonuclear 
kinetics. While detailed discussion of individual rates is beyond the scope of 
this article, the interested reader is directed to the following reviews, in this 
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volume and elsewhere. Experimental nuclear rates and the methodology for 
measuring them have been reviewed in detail by [4,7]. A number of articles in 
this volume review the charged particle reaction rates which are important in 
the Sun [8], in later stages of stellar evolution [9], and for explosive burning 
on compact objects [10,11,12]. Experimental neutron capture cross sections 
are also summarized [13,14]. Theoretical modeling of these rates [15,16] is vi- 
tally important to provide the many rates for which experimental information 
is incomplete or non-existant. Beyond measuring or calculating these rates, 
progress in nuclear astrophysics also requires the compilation and dissemi- 
nation (in a usable form) of this hard-earned nuclear data to the broadest 
astrophysical community [17,18,19]. 

When particle k in Eq. 2 is a photon, the distribution d 3 rik is given by the 
Plank distribution. Furthermore, the relative velocity is always c and thus the 
integral is separable, simplifying to 



T 3 = 



Jd 3 nj J ca(E 7 )E, 



I eJL/LT)-l dE ^ X ^ T)n - (6) 



7r 2 {ch) 3 J o exp(£ 7 /A; B T) - 1 



In practice it is not usually necessary to directly evaluate the photodisinte- 
gration cross sections (see, however, [20] for exceptions), because they can be 
expressed by detailed balance in terms of the capture cross sections for the 
inverse reaction, / + m — » j + 7 [6] . 

A*r(T) = (^)(^) 3/2 ( ! |^) 3/2 (^)exp(-g ta A B T). (7) 



This expression depends on the partition functions, Gk = Z)i(2Jj+l) exp(— Ei/ksT) 
(which account for the populations of the excited states of the nucleus), the 
mass numbers, A, the temperature T, the inverse reaction rate (I, m), and the 
reaction Q- value (the energy released by the reaction), Q im = (mi + m m — 
rrij)c 2 . Since photodisintegrations are endoergic, their rates are vanishingly 
small until sufficient photons exist in the high energy tail of the Planck distri- 
bution with energies > Qi m . As a rule of thumb this requires TttQ lm /30kB- 

In practice, these experimental and theoretical reaction rates are determined 
for bare nuclei, while in astrophysical plasmas, these reactions occur among a 
background of other nuclei and electrons. As a result of this background, the 
reacting nuclei experience a Coulomb repulsion modified from that of bare nu- 
clei. For high densities and/or low temperatures, the effects of this screening of 
reactions becomes very important. Under most conditions (with non- vanishing 
temperatures) the generalized reaction rate integral can be separated into the 
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traditional expression without screening [Eq. 5] and a screening factor, 



This screening factor is dependent on the charge of the involved particles, the 
density, temperature, and the composition of the plasma. For more details on 
the form of f scr , see, e.g., [21,22,23,24]. At high densities and low temperatures 
screening factors can enhance reactions by many orders of magnitude and 
lead to pycnonuclear ignition. In the extreme case of very low temperatures, 
where reactions are only possible via ground state oscillations of the nuclei 
in a Coulomb lattice, Eq. 8 breaks down, because it was derived under the 
assumption of a Boltzmann distribution [25,26]. 

In stellar plasmas, target nuclei also do not exist solely in their ground states. 
This complicates the rate expression in Eq. 5, which now must take into ac- 
count the cross sections for capture out of the different excited states and 
properly weight them according to their probability of occurrence in the en- 
semble of target nuclei. Because the timescales for transitions between excited 
states of a nucleus are typically much shorter than other reaction timescales, 
it is usually valid to assume that the nuclei are internally equilibrated and 
a given excited state is populated in the ensemble by the usual Boltzmann 
factor oc e~ E / fesT , where E now is the excitation energy of that state. From 
this, one may derive a factor, called the stellar enhancement factor (SEF), to 
correct the ground-state reaction rate for the population of excited states (see, 



Interesting complications to this prescription arise when the internal equili- 
bration timescale for a nucleus exceeds the other reaction timescales in the 
problem. This usually occurs when a large spin difference between the ground 
state of an isotope and its first excited state prevents them from communi- 
cating directly via internal transitions. The isotope 26 Al is the prototypical 
example, with a ground and isomeric state of spin and parity 5 + and + , 
respectively. When these two states equilibrate, it is not through a direct 
transition between the states-that simply takes too long. Rather the equi- 
libration occurs through transitions to higher-lying levels that communicate 
effectively with both levels. This problem has been studied in detail (e.g., [29]), 
and the most straightforward way of dealing with this is to assume that the 
higher-lying levels are in a steady state since their destruction timescales are 
so rapid. When this is done, the population of the target nucleus break up into 
two ensembles, one tied to the ground state and one tied to the isomer [30]. 
These may then be treated as separate species in the network, each with its 
own rates to and from other isotopes and with effective rates for transitions 
between the two ensembles. Fortunately, the number of nuclei requiring such 
treatment is small and does not add much computational burden to running 
the nuclear reaction network. 





e.g., [27,28]). 



5 



A procedure similar to that used to derive Eq. 6 applies to captures of elec- 
trons by nuclei. Because the electron is 1836 times less massive than a nucleon, 
the velocity of the nucleus j in the center of mass system is negligible in com- 
parison to the electron velocity (\vj — v e \ ~ \v e \). In the neutral, completely 
ionized plasmas typical of the astrophysical sites of nucleosynthesis, the elec- 
tron number density, n e , is equal to the total density of protons in nuclei, 
Yj% ZiUi. However in many of these astrophysical settings the electrons are at 
least partially degenerate, therefore the electron distribution cannot be as- 
sumed to be Maxwellian. Instead the capture cross section may be integrated 
over a Boltzmann, partially degenerate, or degenerate Fermi distribution of 
electrons, depending on the astrophysical conditions. The resulting electron 
capture rates are functions of T and n e , rj = Xj >e (T,n e )nj. Similar equations 
apply for the capture of positrons which are in thermal equilibrium with pho- 
tons, electrons, and nuclei. Electron and positron capture calculations have 
been performed for a large variety of nuclei with mass numbers between A=20 
and A=100 (see [31] for more information). For normal decays, like beta or 
alpha decays, with a characteristic half-life T1/2, Eq. 6 also applies, with the 
decay constant Xj = ln2/ri/ 2 . In addition to innumerable experimental half- 
life determinations, beta-decay half-lives for a wide range of unstable nuclei 
have been predicted (see [32,33]). 

Even though the size of the neutrino scattering cross section on nuclei and 
electrons is very small, at high densities (p ~ 10 13 gcm -3 ), enough scattering 
events occur to thermalize the neutrino distribution. Under such conditions 
the inverse process to electron capture (neutrino capture) can occur in sig- 
nificant numbers and the neutrino capture rate can be expressed in a form 
similar to Eqs. 6 by integrating over the thermal neutrino distribution (e.g. 
[34]). Inelastic neutrino scattering on nuclei can also be expressed in this form. 
The latter can cause particle emission, similar to photodisintegration. In this 
volume, see [35,36,37] for more complete discussion of the interactions of neu- 
trinos with matter. The calculation of these rates can be further complicated 
by the neutrinos being out of thermal equilibrium with the local environment. 
When thermal equilibrium among neutrinos was established at a different lo- 
cation, then the neutrino distribution might be characterized by a chemical 
potential and a temperature different from the local values. Otherwise, the 
neutrino distribution must be evolved in detail (see, e.g., [38]). 



2.2 Thermonuclear Rate Equations 

The large number of reaction types discussed in §2.1 can be divided into 
3 functional categories based on the number of reactants which are nuclei. 
The reactions involving a single nucleus, which include decays, electron and 
positron captures, photodisintegrations, and neutrino induced reactions, de- 
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pend on the number density of only the target species. For reaction involving 
two nuclei, the reaction rate depends on the number densities of both target 
and projectile nuclei. There are also a few important three-particle process 
(like the triple-a process, see §5) which are commonly successive captures 
with an intermediate unstable target (see, e.g., [39,40]). Using an equilibrium 
abundance for the unstable intermediate, the contributions of these reactions 
are commonly written in the form of a three-particle processes, depending 
on a trio of number densities. Grouping reactions by these 3 functional cate- 
gories, the time derivatives of the number densities of each nuclear species in 
an astrophysical plasma can be written in terms of the reaction rates, r, as 



= »>, • • I.\;,,r,,. (9) 

p=const j j,k j,k,l 



where the three sums are over reactions which produce or destroy a nucleus 
of species % with 1, 2 & 3 reactant nuclei, respectively. The M s provide for 
proper accounting of numbers of nuclei and are given by: Af* = N iy Af^ k = 

Ayrim=il^m|!, and A/£ fc) , = ^/rim=i \N m \l The N!s can be positive or 
negative numbers that specify how many particles of species i are created or 
destroyed in a reaction, while the denominators, including factorials, run over 
the n,j t k or rij t k,i different species destroyed in the reaction and avoid double 
counting of the number of reactions when identical particles react with each 
other (for example in the 12 C + 12 C or the triple-a reactions; for details see 
[6])- 

In addition to nuclear reactions, expansion or contraction of the plasma can 
also produce changes in the number densities rij. To separate the nuclear 
changes in composition from these hydrodynamic effects, the nuclear abun- 
dance Yi = rii/pN A , where N A is Avogadro's number, is commonly used. For 
a nucleus with atomic weight A^Yi represents the mass fraction of this 
nucleus, therefore = 1. Likewise, the equation of charge conservation 

becomes Ztfi — Y e , where Y e (— n e / pN A ) is the electron abundance (also 
termed the electron fraction). By recasting Eq. 9 in terms of nuclear abun- 
dances Yi, a set of ordinary differential equations for the evolution of Yj results 
which depends only on nuclear reactions. In terms of the reaction cross sec- 
tions introduced in §2.1, this reaction network is described by the following 
set of differential equations 

Yi = + J2Afl kP N A (j, k)Y 3 Y k + J2Af; Al p 2 N 2 A (j, k, o w^(io) 

3 j,k j,k,l 
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2.3 Coupling Nuclear Reaction Networks to Hydrodynamics 

As we touched on in the previous section, nuclear processes are tightly linked 
to the hydrodynamic behavior of the bulk medium. Thermonuclear processes 
release (or absorb) energy, which alters the pressure and causes hydrodynamic 
motions. These motions may disperse the thermonuclear ash and bring a con- 
tinued supply of fuel to support the thermonuclear flame. The compositional 
changes caused by thermonuclear reactions can also change the equation of 
state and opacity, further impacting the hydrodynamic behavior. For purposes 
of this review of thermonuclear kinetics methods, which generally assume that 
thermonuclear and hydrodynamic changes in local composition can be success- 
fully decoupled (or treated in an operator split fashion) , we include only a brief 
description of how this decoupling is best achieved. Miiller [41] provides an 
authoritative overview and discusses the difficulties (and open issues) involved 
when including nucleosynthesis within hydrodynamic simulations. 

The coupling between thermonuclear processes and hydrodynamic changes can 
be divided into two categories by considering the spatial extent of the cou- 
pling. Nucleosynthetic changes in composition and the resultant energy release 
produce local changes in hydrodynamic quantities like pressure and temper- 
ature. The strongest of these local couplings is the release (or absorption) of 
energy and the resultant change in temperature. Changes in temperature are 
particularly important because of the exponential nature of the temperature 
dependence of thermonuclear reaction rates. Since the nuclear energy release 
is uniquely determined by the abundance changes, the rate of thermonuclear 
energy release, e, is given by 

e nuc = ~ E N A M lC 2 Y t ( MeVg- 1 s" 1 ). (11) 



where M^c 2 is the rest mass energy of species % in MeV. Since all reactions con- 
serve nucleon number, the atomic mass excess M ex ^ = Mj — Aiin u (m u is the 
atomic mass unit) can be used in place of the mass Mj in Eq. 11 (see [42] for 
a recent compilation of mass excesses). The use of atomic mass units has the 
added benefit that electron conservation is correctly accounted for in the case 
of f3~ decays and e~ captures, though reactions involving positrons require 
careful treatment. In general, the nuclear energy release is deposited locally, 
so the rate of thermonuclear energy release is equal to the nuclear portion of 
the hydrodynamic heating rate. However, there are instances where nuclear 
products do not deposit their energy locally. Escaping neutrinos can carry 
away a portion of the thermonuclear energy release. In the rarefied environ- 
ment of supernova ejecta at late times, positrons and gamma rays released by 
(3 decays are not completely trapped. In most such cases, the escaping particles 
stream freely from the reaction site, allowing adoption of a simple loss term 
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analogous to Eq. 11 with MjC 2 replaced by an averaged energy loss term. For 
this reason, weak reaction rate tabulations provide averaged neutrino losses. 
From these we can construct 

loss ^ ^ (-^v) i,weaki (12) 



where we consider only those contributions to Y due to neutrino producing 
reactions. In some cases, like supernovae, subsequent interactions between the 
escaping leptons or gamma rays and matter require more complete transport 
to be considered. Other important quantities which are impacted by nucle- 
osynthesis, like Y e , can be obtained by appropriate sums over the abundances 
and also need not be evolved separately. 

Implicit solution methods require the calculation of Y(t + At), where At is 
the nuclear timestep, which in turn requires knowledge of T(t + At). One 
could write a differential equation for the energy release analogous to Eq. 10, 
with the J\fs replaced by the reaction Q-values, and thereby evolve the energy 
release (and calculate temperature changes) as an additional equation within 
the network solution. Miiller [44] has shown that such a scheme can help avoid 
instabilities in the case of a physically isolated zone entering or leaving nuclear 
statistical equilibrium. In general, however, use of this additional equation 
is made unnecessary by the relative slowness with which the temperature 
changes. The timescale on which the temperature changes is given by 

r T = T/f » C v T/e nuc (13) 



and is often called the ignition timescale. The timescale on which an individual 
abundance changes is its burning time, 

r( A Z) = F( A Z)/F( A Z) ~ minr fe ( A Z) (14) 



where r^( A Z) is defined in Eq. 4. In general tt differs from t( a Z) of the 
principle fuel by the ratio of thermal energy content to the energy released 
by the reaction. For degenerate matter this ratio can approach zero, allowing 
for explosive burning. In contrast, accurate prediction of less abundant, but 
still important, species requires that the reaction network timestep At be 
chosen to be the burning timescale of a less abundant species, typically with 
an abundance of 10~ 6 or smaller [45]. Since the dominant fuel is most often one 
of the more abundant constituents and the burning timescales are proportional 
to the abundance, tt is typically an order of magnitude or more larger than 
the reaction network timestep (see, e.g., [46,47]). It is therefore sufficient to 
calculate the energy gain at the end of a timestep via equation 11, modified 
as discussed above, and approximate T(t + St) ps T(t) or to extrapolate based 
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on e(t). Since other locally coupled quantities have characteristic timescales 
much longer than At, they too can be decoupled in a similar fashion. For 
the remainder of this review, we will consider only the equations governing 
changes in isotopic abundances, remembering that additional equations can 
easily be constructed for those special circumstances where they are necessary. 

Spatial coupling, particularly the modification of the composition by hydrody- 
namic movements such as diffusion, convective mixing and advection (in the 
case of Eulerian hydrodynamics methods), represents a more difficult chal- 
lenge. By necessity, an individual nucleosynthesis calculation examines the 
abundance changes in a locality of uniform composition. The difficulties as- 
sociated with strong spatial coupling of the composition occur because this 
nucleosynthetic calculation is spread over an entire hydrodynamic zone. Con- 
vection can result in strong abundance gradients across a single hydrodynamic 
zone, which with the assumption of compositional uniformity, can result in 
very different outcomes as a function of the fineness of the hydrodynamic 
grid. Thermonuclear supernovae, where the thickness of the flame front can 
be billionths of the radius of the white dwarf [48,49], present an extreme exam- 
ple of microscopic nuclear inhomogeneity. Eulerian advection of compositional 
boundaries can also have extremely unphysical consequences. Fryxell et al. [50] 
demonstrated how this artificial mixing can produce an unphysical detonation 
in a shock tube calculation by mixing cold unburnt fuel into the hot burnt re- 
gion. A related problem is the conservation of species. Hydrodynamic schemes 
must carefully conserve the abundances (or partial densities) of all species 
[50,51,52], lest they provide unphysical abundances to the nucleosynthesis cal- 
culations, which must assume conservation, and thereby produce unphysical 
results. Because of these problems, hydrodynamic methods with excellent cap- 
ture of shocks and contact or compositional discontinuities are best suited to 
nucleosynthesis calculations. 

The relative size of the burning timescales, when compared to the relevant 
diffusion, sound crossing or convective timescale, dictates how these problems 
must be addressed. If all of the burning timescales are much shorter than 
the timescale on which the hydrodynamics changes the composition, then 
the assumption of uniform composition is satisfied and the nucleosynthesis of 
each hydrodynamic zone can be treated independently. If all of the burning 
timescales are much larger than, for example, the convective timescale, then 
the composition of the entire convective zone can be treated as uniform and 
slowly evolving. The greatest complexity occurs when the timescales on which 
the hydrodynamics and nucleosynthesis change the composition are similar. 
Oxygen shell burning represents an excellent example of this as the sound 
travel, convective turnover and nuclear burning timescales are all of the same 
order as the evolutionary time. The results of 2D simulations [53] demonstrate 
convective overshooting, highly non-uniform burning and a velocity structure 
dominated by convective plumes. In section 3.1 we will briefly discuss the nu- 
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merical challenges posed by fully coupling nucleosynthesis and nuclear mixing. 
Silicon burning represents a different challenge [5], as the timescales for the 
transformation of silicon to iron are much slower than the convective turnover 
time, but the burning timescales for the free neutrons, protons and a-particles 
which maintain QSE are much faster, providing a strong motivation for the 
hybrid reaction networks we will discuss in §5. 



3 Solving the Nuclear Reaction Network 

In principle, the initial value problem presented by the nuclear reaction net- 
work can be solved by any of a large number of methods discussed in the 
literature. However the physical nature of the problem, reflected in the A's 
and (cn>)'s, greatly restricts the optimal choice. The large number of reactions 
display a wide range of reaction timescales, r (see Eq. 4). Numerical systems 
whose solutions depend on a wide range of timescales are termed stiff. Gear 
[54] demonstrated that even a single equation can be stiff if it has both rapidly 
and slowly varying components. Practically, stiffness occurs when the limita- 
tion of the timestep size is due to numerical stability rather than accuracy. 

— * 

— * 

A more rigorous definition [55] is that a system of equations Y(Y) is stiff if 

— * 

— # 

the eigenvalues A-,- of the Jacobian dY/dY obey the criterion that for negative 
9ft(Aj) (the real part of the eigenvalues \j) 



for j — 1, • • • , N. As we will explain in this section, S > 10 15 is not uncommon 
in astrophysics. 

Astrophysical calculations of nucleosynthesis belong to the general field of re- 
active flows, and therefore share some characteristics with related terrestrial 
fields. In particular, chemical kinetics, the study of the evolution of chemical 
abundances, is an important part of atmospheric and combustion physics and 
produces sets of equations much like Eq. 9 (see [56] for a through introduction). 
These chemical kinetics systems are known for their stiffness and a great deal 
of effort has been expended on developing methods to solve these equations. 
Many of the considerations for the choice of solution method for chemical 
kinetics also apply to thermonuclear kinetics. In both cases, temporal integra- 
tion of the reaction rate equations is broken up into short intervals because 
of the need to update the hydrodynamics variables. This favors one step, self 
starting algorithms. Because abundances must be tracked for a large number 
of computational cells (hundreds to thousands for one dimensional models, 
millions to billions for the coming generation of three dimensional models), 




(15) 
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memory storage concerns favor low order methods since they don't require the 
storage of as much data from prior steps. In any event, both the errors in fluid 
dynamics and in the reaction rates are typically a few percent or more, so the 
greater precision of these higher order methods often does not result in greater 
accuracy. Note that these statements do not in general apply to calculations 
of Big Bang nucleosynthesis (at least those that assume homogeniety) . As a 
result quite different methods are employed in these calculations [57] 

Because of the wide range in timescales between strong, electromagnetic and 
weak reactions, nuclear reaction networks are extraordinarily stiff. PP chain 
nucleosynthesis, responsible for the energy output of the Sun, offers an ex- 
cellent example of the difficulties. The first reaction of the PP1 chain is 
1 H(p, e + z/) 2 H, the fusion of two protons to form deuterium. This is a weak 
reaction, requiring the conversion of a proton into a neutron, and releasing a 
positron and a neutrino. As a result, the reaction timescale r p ( 1 H) is very long, 
billions of years for conditions like those in the solar interior. The second reac- 
tion of the PP1 chain is the capture of a proton on the newly formed deuteron, 
2 H(p, 7) 3 He. For conditions like those in the solar interior, the characteristic 
timescale, r p ( 2 H) is a few seconds. Thus the timescales for two of the most 
important reactions for hydrogen burning in stars like our Sun differ by more 
than 17 orders of magnitude (see [3,8] for a more complete discussion of the PP 
chain). This disparity results not from a lack of 1 H + 1 H collisions (which oc- 
cur at a rate F( 1 H)/F( 2 H) ~ 10 17 times more often than 1 H + 2 H collisions), 
but from the rarity of the transformation of a proton to a neutron. While the 
presence of weak reactions among the dominant energy producing reactions is 
unique to hydrogen burning, most nucleosynthesis calculations are similarly 
stiff, in part because of the need to include weak interactions but also the 
potential for neutron capture reactions, which occur very rapidly even at low 
temperature, following any release of free neutrons. The nature of the nuclear 
reaction network equations has thus far limited the astrophysical usefulness of 
the most sophisticated methods to solve stiff equations developed for chem- 
ical kinetics. However work to harness this resource continues. For example, 
Timmes [58] has examined higher order Kaps-Rentrop and Bader-Deuflhard 
semi-implicit time integration algorithms while Mott et al. [59] have studied 
asymptotic and quasi-steady state methods, which do not require a matrix 
solution. 

For a set of nuclear abundances Y, one can calculate the time derivatives of 

the abundances, Y using Eq. 10. The desired solution is the abundance at a 
future time, Y(t + At), where At is the network timestep. For simplicity, most 
past and present nucleosynthesis calculations use the simple finite difference 
prescription 

Y{t + At) ~ Y (t) = (1 - Q)Y(t + At) + QY(t). (16) 
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This choice is also supported by the advantages of coupling low order, single 
step methods with hydrodynamics. With = 1, Eq. 16 becomes the explicit 
Euler method while for O = it is the implicit backward Euler method, both 
of which are first order accurate. For O = 1/2, Eq. 16 is the semi-implicit 
trapezoidal method, which is second order accurate. For the stiff set of non- 
linear differential equations which form most nuclear reaction networks, a fully 
implicit treatment is generally most successful [45], though the trapezoidal 
method has been used in Big Bang nucleosynthesis calculations [60], where 
coupling to hydrodynamics is less important. Solving the fully implicit version 
of Eq. 16 is equivalent to finding the zeros of the set of equations 

g(t + 5 + ^ + = 0. (17) 



This is done using the Newton-Raphson method (see, e.g., [61]), which is 
based on the Taylor series expansion of Z(t + At), with the trial change in 
abundances given by 

AY=M t + At) Yz, (18) 
\dY(t + At)J 1 ; 



where dZj dY is the Jacobian of Z. 

Historically [45] and in some modern applications (see, e.g., [62]), each timestep 
consists of only a single application of the procedure outlined in Eqs. 17 & 18. 
This semi-implicit backward Euler method has the advantage of a relatively 
small and predictable number of matrix solutions, but there are only heuristic 
checks that the chosen timestep results in a stable or accurate solution. In fully 
implicit backward Euler schemes, iteration using the procedure of Eqs. 17 & 

— * — * 

18 continues until both AY and Z are below some tolerance, providing a mea- 
sure of the stability and accuracy. If this convergence does not occur within 
a reasonable number of iterations, the timestep is subdivided into smaller in- 
tervals until a converged solution can be achieved, allowing the fully implicit 
backward Euler integration to respond to instability or inaccuracy in a way 
that is impossible with the semi-implicit backward Euler approach. Higher 
order methods allow better estimates of the truncation error by comparing 
the solutions of different order schemes and sub-dividing the timestep if these 
errors are too large. 

A potential numerical problem with the solution of Eq. 17 is the singularity 
of the Jacobian matrix, dZ(t + At)/dY(t + At). From Eq. 17, the individual 
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matrix elements of the Jacobian have the form 



dZi Sij ())', Sij 




(19) 



dYj At dYj At 



where 5^ is the Kronecker delta, and Tj(i) is the destruction timescale of nu- 
cleus i with respect to nucleus j for a given reaction, as defined in Eq. 4. The 
sum accounts for the fact that there may be more than one reaction by which 
nucleus j is involved in the creation or destruction of nucleus %. Along the 
diagonal of the Jacobian, there are two competing terms, 1/ At and X) 1/Tj(i). 
This sum is over all reactions which destroy nucleus i, and is dominated by the 
fastest reactions. As a result, J2 l/ r i(*) can be orders of magnitude larger than 
the reciprocal of the desired timestep, 1/At. This is especially a problem near 
equilibrium, where both destruction and the balancing production timescales 
are very short in comparison to the preferred timestep size, resulting in differ- 
ences close to the numerical accuracy (i.e. 14 or more orders of magnitude). In 
such cases, the term 1/ At is numerically neglected, leading to numerically sin- 
gular matrices. One approach to avoiding this problem is to artificially scale 
these short, equilibrium timescales by a factor which brings their timescale 
closer to At, but leaves them small enough to ensure equilibrium. While this 
approach has been used successfully, the ad hoc nature of this artificial scaling 
renders these methods fragile. A more promising approach is to make directly 
use of equilibrium expressions for abundances, which, as we will discuss in §5, 
also ensures the economical use of computer resources. 

3.1 Taking Advantage of Matrix Sparseness 

For larger nuclear reaction networks, the Newton-Raphson method requires 
solution of a moderately large (N = 100 — 3000) matrix equation for each 
zone. Since general solution of a dense matrix scales as 0(N 3 ), this can make 
these large networks progressively much more expensive. While in principal, 
every species reacts with each of the hundreds of others, resulting in a dense 
Jacobian matrix, in practice it is possible to neglect most of these reactions. 
Because of the ZiZj dependence of the repulsive Coulomb term in the nuclear 
potential, captures of free neutrons and isotopes of H and He on heavy nu- 
clei occur much faster than fusions of heavier nuclei. Furthermore, with the 
exception of the Big Bang nucleosynthesis and PP-chains, reactions involv- 
ing secondary isotopes of H (deuterium and tritium) and He are neglectable. 
Likewise, photodisintegrations tend to eject free nucleons or a-particles. Thus, 
with a few important exceptions, for each nucleus we need only consider twelve 
reactions linking it to its nuclear neighbors by the capture of an n, p, a or 7 
and release a different one of these four. The exceptions to this rule are the 
few heavy ion reactions important for burning stages like carbon and oxy- 
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Fig. 1. Graphic demonstration of the sparseness of the Jacobian matrix. The filled 
squares represent the non-zero elements. 

gen burning where the dearth of light nuclei cause the heavy ion collisions to 
dominate. 

Fig. 1 demonstrates the sparseness of the resulting Jacobian matrix, for a 300 
nuclei network chosen to handle all the energy generating stages in the life 
of a massive star. The nuclei are indexed in order of increasing Z and then 
A. Of the 90,000 matrix elements, fewer than 5,000 are non-zero. In terms 
of the standard forms for sparse matrices, this Jacobian is best described as 
doubly bordered, band diagonal. With a border width, AB, of 45 necessary 
to include the heavy ion reactions among 12 C, 16 and 20 Ne along with the 
free neutrons, protons and a-particles and a band diagonal width, AD, of 
54, even this sparse form includes almost 50,000 elements. With solution of 
the matrix equation consuming 90+% of the computational time, there is 
clearly a need for custom tailored solvers which take better advantage of the 
sparseness of the Jacobian [43]. To date best results for small (N<100) matrices 
are obtained with machine optimized dense solvers (e.g. LAPACK) or matrix 
specific solvers generated by symbolic processing [41,44]. For large matrices, 
generalized sparse solvers, both custom built and from software libraries, are 
used (see, [58]). 

Thus far in this section we have considered thermonuclear kinetics in the limit 
where the hydrodynamics is treated separately by operator-splitting. How- 
ever, as we discussed in §2.3, this is approach is not always viable, so we 
conclude this section with some remarks on the additional challenges encoun- 
tered in recent efforts to solve coupled nucleosynthesis and hydrodynamical 
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mixing. Because these calculations must also be finite differenced with an im- 
plicit scheme, they too require matrix solutions and all zones must be solved 
concurrently. The resulting matrices are quite large: without consideration of 
matrix sparseness, one would need to store N 2 M 2 numbers, where N is the 
number of species in each zone and M is the number of zones. Fortunately, 
however, the resulting matrices are quite sparse since they are block diagonal 
(each block is similar in form to that shown in Fig. 1) with bands linking adja- 
cent spatial zones. Such large sparse matrices are typically solved by iterative 
techniques, such as biconjugate gradient which requires only matrix multiplies 
and, thus, only storage of the non-zero matrix elements. To date, computation 
of iV = 1000, M = 1000 problems (that is, 10 6 x 10 6 matrices) can be fairly 
routinely handled on a single processor. This approach to solving coupled nu- 
cleosynthesis and nuclear mixing is necessary in cases where operator splitting 
is suspect, and we expect to see it finding greater use in future calculations. 

3.2 Physically Motivated Network Specialization 

Often from a physical understanding one can specialize the general solution 
method and thereby greatly reduce the computational cost. As an example, 
we briefly discuss the r-process approximation described in [63,64]. For nuclei 
with A > 100, charged particle captures (proton and a) as well as their reverse 
photodisintegrations virtually cease when T < 3 GK. This leaves only neutron 
captures and their reverse photodisintegration reactions, as well as /9-decays, 
which can also lead to the emission of delayed neutrons. In this case, Eq. 10 
greatly simplifies, leaving 

r( A z) = ^M^y^z) + a^ +1 f( a+1 z) + f Ag? M+i y( A +Jz - 1) 

3=0 

~ (n n (o~v) n £ A + \\ A + £ A§5 j Y ( A Z) , (20) 

where (cv) n z\ and X% A are the velocity integrated neutron capture cross sec- 
tion and the photodisintegration rate for the nucleus A Z, while X^za) * s ^ e 
cay constant for the f3~ decay of A Z, with j delayed neutrons (up to a maximum 
of J). The assumption is made that the neutron abundance (Y n = n n /pNA) 
varies slowly enough that it may be evolved explicitly. One can see that in 
Eq. 20, with n n thereby fixed, the time derivatives of each species have a lin- 
ear dependence on only the abundances of their neighbors in the same isotopic 
chain (nuclei with the same Z), or that with one less proton (Z — 1). One can 
then divide the network into separate pieces for each isotopic chain, and solve 
them sequentially, beginning with the lowest Z. The "boundary" terms for this 
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lowest Z chain can be supplied by a previously run or concurrently running 
full network calculation which need extend only to this Z . This reduces the 
solution of a matrix with more than a thousand rows to the solution of roughly 
30 smaller matrices. Furthermore each of these smaller matrices is also tridi- 
agonal increasing speed further. Freiburghaus et al. [64] tested the assumption 
of slow variation in the neutron abundance and have demonstrated the use- 
fulness of this method in r-process simulations, achieving a large decrease in 
computational cost. A similar treatment has been successfully applied to ex- 
plosive hydrogen burning based on the assumption of slowly varying proton 
and alpha abundances [65]. As we will discuss in §5, for other burning stages 
there also exist physically motivated simplifications to the general network 
solution method. 



4 Equilibria in Nuclear Astrophysics 

An isolated thermodynamic system tends to evolve towards an equilibrium 
state, and whether the system reaches that equilibrium depends on the compe- 
tition between the timescale for equilibration and the other relevant timescales. 
For example, the reaction-rate expression found in Eq. (5) assumed a local 
thermal equilibrium in the plasma such that the velocities of the interacting 
particles are well described by a Maxwell-Boltzmann distribution, which is de- 
termined by a single parameter, the local temperature. This is valid because 
the timescale to achieve this equilibrium is much shorter than the other inter- 
action timescales in the plasma. Similarly, as discussed in §2.1, it is usually 
assumed in nuclear astrophysics that the internal states of a nucleus are also 
in equilibrium so that the probability of finding a nucleus in particular excited 
state is proportional to the Boltzmann factor for that state, which, as before, is 
determined by the temperature and the (presumably known) excitation energy 
of the state. Again this is generally a valid assumption because the transitions 
within the nucleus are typically much faster than other interactions, although 
interesting challenges arise when this assumption is not valid (e.g., [29,30]). 

At conditions of high temperature and density, the thermonuclear reaction 
rates themselves may be sufficiently rapid to achieve equilibrium within the 
timescale set by the hydrodynamics of the astrophysical setting. This permits 
considerable simplification of the calculation of the nuclear abundances. In 
most such cases, the fast strong and electromagnetic reactions reach equilib- 
rium while those involving the weak nuclear force do not. Since the weak reac- 
tions are not equilibrated, the resulting Nuclear Statistical Equilibrium (NSE) 
requires monitoring of weak reaction activity. Even with this stricture, NSE 
offers many advantages, since hundreds of abundances are uniquely defined 
by the thermodynamic conditions and a single measure of the weak interac- 
tion history or the degree of neutronization. Computationally, this reduction 
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in the number of independent variables greatly reduces the cost of nuclear 
abundance evolution. Because there are fewer variables to follow within a hy- 
drodynamic model, the memory footprint of the nuclear abundances is also 
reduced, an issue of importance in modern mult i- dimensional models. Finally, 
the equilibrium abundance calculations depend on binding energies and par- 
tition functions, quantities which are better known than many reaction rates. 
This is particularly true for unstable nuclei and for conditions where the mass 
density approaches that of the nucleus itself, resulting in exotic nuclear struc- 
tures. 

The expression for NSE is commonly derived using either chemical potentials 
or detailed balance (see, e.g., [3]). For a nucleus A Z, composed of Z protons and 
N = (A — Z) neutrons, in equilibrium with these free nucleons, the chemical 
potential of A Z can be expressed in terms of the chemical potentials of the 
free nucleons 

fi z ,A = Zfx p + Np n . (21) 



Substituting the expression for the Boltzmann chemical potential (including 
rest mass) into Eq. 21 allows derivation of an expression for the abundance of 
every nuclear species in terms of the abundances of the free protons (Y p ) and 
neutrons (Y n ), 



where G{ A Z) and B( A Z) are the partition function and binding energy of the 
nucleus A Z, Na is Avagadro's number, ks is Boltzmann's constant, p and T 
are the density and temperature of the plasma, and 6 = (m u k B T /2-irh 2 ) 3 / 2 . 

Abundances of all nuclear species can therefore be expressed as functions of 
two quantities. Nucleon number conservation (J2 AY = 1) provides one con- 
straint. The second constraint is the amount of weak reaction activity, often 
expressed in terms of the total proton abundance, X ZY, which charge con- 
servation requires equal the electron abundance, Y e . Thus the nuclear abun- 
dances are uniquely determined for a given (T,p,Y e ). Alternately, the weak 
interaction history is sometimes expressed in terms of the neutron excess 
7] = XX iV — Z)Y . Figure 2 displays the temperature and density dependence of 
A = X AY I J2 Y = 1/ J2 Y, the average nuclear mass of the NSE distribution. 
At high temperatures, free nucleons are favored, hence A ~ 1. For intermediate 
temperatures the compromise of retaining large numbers of particles while in- 
creasing binding energy favors 4 He, which has 80% of the binding energy of the 
iron peak nuclei. At low temperatures, Eq. 22 strongly favors the most bound 
nuclei, the iron peak nuclei, so A — > 60 as the temperature drops. Density can 
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Fig. 2. The average atomic mass for material in NSE as a function of Temperature. 
The solid lines include screening corrections to the nuclear binding energies, while 
the dotted lines ignore this effect. 



be seen to scale the placement of these divisions between high, intermediate 
and low temperature and increase the average mass at low temperature. Vari- 
ations in Y e do not strongly affect Figure 2. At high temperatures, it simply 
effects the ratio of Y p /Y n m Y e /1 — Y e . At low temperatures, variation in Y e 
changes which Fe-peak isotopes dominate. For example, though 56 Ni is less 
tightly bound than 54 Fe, it is more tightly bound than 54 Fe + 2 1 H, which 
would be required by charge conservation if Y e ~ .5. Thus F( 56 Ni) > F( 54 Fe) 
for low T with Y e ~ .50, but F( 54 Fe) > F( 56 Ni) for smaller Y e . In general, 
the most abundant nuclei at low temperatures are the most bound nuclei for 
which Z/A ~ Y e . 

As with any equilibrium distribution, there are limitations on the applicability 
of NSE. For NSE to provide a good estimate of the nuclear abundances the 
temperature must be sufficient for the endoergic reaction of each reaction pair 
to occur. Since for all particle-stable nuclei between the proton and neutron 
drip lines (with the exception of nuclei unstable against alpha decay), the 
photodisintegrations are endoergic, with typical Q-values among (f3) stable 
nuclei of 8-12 MeV, by Eq. 7 this requirement reduces to T > 3GK. While 
this requirement is necessary, it is not sufficient. In the case of hydrostatic 
silicon burning, even when this condition is met, appreciable time is required 
to convert Si to Fe-peak elements. In the case of explosive silicon burning, 
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the adiabatic cooling on timescales of seconds can cause conditions to change 
more rapidly than NSE can follow, breaking down NSE first between 4 He and 
12 C, at T ~ 6GK [66] and later between the species near silicon and the 
Fe-peak nuclei, at T ~ 4GK [67,68]. Thus it is clear that in the face of suffi- 
ciently rapid thermodynamic variations, NSE provides a problematic estimate 
of abundances. It is important to note that, in spite of the breakdown of the 
global NSE, many nuclei under these conditions do obey local equilibrium 
relations with their neighbors. 



5 Merging Equilibria with Reaction Networks 

In spite of the limitations on the applicability of NSE, the reduced compu- 
tational cost provides a strong motivation to maximize the use of equilibria. 
The use of equilibrium expressions for single abundances is, in fact, common in 
nuclear reaction networks, typically to track the abundances of short-lived un- 
stable intermediates in "three-particle" processes. The most common example 
of this is the triple a process, 

4 He + 4 He ^ 8 Be Q = -.09 MeV 

(23) 

4 He + 8 Be ^ 12 C* -> 12 C + 7 Q = 7.37 MeV , 

by which Helium burning occurs. With r( 8 Be) ~ 10~ 16 sec, only rarely does 
a 8 Be survive long enough for a second a to capture. As a result of the near 
balance of the first reaction pair, the abundance of 8 Be can be expressed in 
terms of the a-particle abundance, 

M(8 y M ° ) VI (24) 



Likewise for temperatures in excess of . 1 GK, the most likely result following 
the second a capture to form an excited state of 12 C is a decay back to 8 Be 
(r a ( 12 C*)/r 7 ( 12 C*) > 10 3 ), thus the abundance of 12 C* is well characterized 
by 

When this is the case, the effective triple a reaction rate is simply that of the 
decay of 12 C from the excited state to the ground state, 

r 3a = pN A YC 2 C*)T^ l2 C*)/h . (26) 



^>^Gfex P 
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This use of local equilibrium within a rate equation shares many characteristics 
with the more elaborate schemes we will discuss next. The number of species 
tracked by the network is reduced since Y( 8 Be) need not be directly evolved. 
Problematically small timescales like r( 8 Be) are removed, replaced by larger 
time scales (t 3q , ~ 10 5 — 10 7 years during core helium burning). The non- 
linearity of network time derivatives is increased (Y( 12 C) oc Y£) under this 
scheme. This approximation also breaks down at low T (for details see [39]). 



6 The QSE-reduced Network 



As we noted in the previous section, while global NSE may not always apply for 
temperatures in the range 3-6 GK, many nuclei are in local equilibrium with 
their neighbors. Beginning with Bodansky et al. [69], a number of attempts 
have been made to take advantage of these partial equilibria (termed quasi- 
equilibria or QSE) to reduce the number of independent variables evolved via 
rate equations and thereby reduce the computational cost of modeling these 
burning stages. To evolve the abundances of every member of a QSE group, 
it is sufficient to evolve the abundance of any single group member along 
with the abundances of the free nucleons. One can thereby reduce the number 
of abundances that are evolved, while still calculating, from QSE relations, 
the abundances of all members of a QSE group and the resultant reaction 
rates, including the electron and neutrino capture reactions responsible for 
changes in the neutronization. The result is a more computationally efficient 
method that retains the accuracy of the full network and yields abundances 
for all nuclei found in the full network. Such methods have been applied to 
the a-network and produced networks that are twice as fast as the minimal 
a-chain network, without significantly affecting the nuclear evolution [70,71]. 
Reductions of an order of magnitude in computational cost have been achieved 
[72,73] for QSE-reduced networks of the size necessary to capture the essential 
features of supernova nucleosynthesis and we believe greater savings are possi- 
ble with further refinement. In addition to silicon burning, there are a number 
of astrophysically important situations where T > Q/30k B for at least some 
of the relevant reactions and so large equilibrium groups exist, but NSE is 
not globally valid. This includes the r-process [63] and the rp-process in X-ray 
bursts (and possibly novae) [65], where neutron or proton separation energies 
(Q n or Qp) of 2 MeVand less are often encountered. 

As a example, we will discuss the QSE-reduced a-network. Its mission is to 
evolve the abundances of the full 14 elements of a conventional a-network 
(which we'll call Y^), and calculate the resulting energy generation, in a more 
efficient way. Under conditions where QSE applies, the existence of the silicon 
and iron peak QSE groups (which are separated by the nuclear shell closures 
Z=N=20 and the resulting small Q-values and reaction rates) allows calcu- 
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lation of these 14 abundances from 7. For the members of the silicon group 
( 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti) and the iron peak group ( 48 Cr, 52 Fe, 56 Ni, 60 Zn) 
the individual abundances can be calculated by expressions similar to Eq. 22, 



Y Q seM A Z) = ^f[)F( 56 Ni)l^, (27) 

where C( Z) is defined in Eq. 22 and {A— 28)/4 and (A— 56)/ 4 are the number 
of a-particles needed to construct A Z from 28 Si and 56 Ni, respectively. Where 
QSE applies, Y T is a function of the abundances of a reduced nuclear set, 71, 
defined as a, 12 C, 16 0, 20 Ne, 24 Mg, 28 Si, 56 Ni and we need only evolve Y n . It 
should be noted that 24 Mg is ordinarily a member of the silicon QSE group 
[5,67,74], but for easier integration of prior burning stages with a conventional 
nuclear network, we will evolve 24 Mg independently. The main task when 
applying such hybrid schemes is finding the boundaries of QSE groups and 
where individual nuclei have to be used instead. Treating marginal group 
members as part of a group increases the efficiency of the calculation, but 
may decease the accuracy. 

While Y n is a convenient set of abundances for calculating Y T ', it is not the 
most efficient set to evolve, primarily because of the non-linear dependence on 
Y a . Instead we define a set of group abundances, Y g = [Y a c, Y( 12 C), F( 16 0), 
F( 20 Ne), F( 24 Mg), Y StG , Y FeG ] where 

idSi group " i£Fe group 

YsiG — E Y t) (28) 

i&Si group 

Y FeG = E Y { . 

iGFe group 

Physically, Y a c represents the sum of the abundances of free a-particles and 
those a-particles required to build the members of the QSE groups from 28 Si 
or 56 Ni, while Y$iG and Yp eG represent the total abundances of the silicon 
and iron peak QSE groups. This method, which here is applied only to the 
chain of a-nuclei can also be generalized to arbitrary networks [73] . For larger 
networks which contain nuclei with N ^ Z, one must be able to follow the 
abundances of free neutrons and protons, particularly since weak interactions 
will change the global ratio of neutrons to protons. In place of Y aG in Eq. 28, 
one constructs Y NG = Y,i,u g ht N iYi + T,i,si( N i ~ 14 )^ + T,i,Fe( N i ~ 2 %) Y i and 
Yzg = flight Z*Y, t + Z itSi (Zi - U)Yi + Ei, Fe ( Z i ~ 28 ) F - if 2§Si and 56Ni are 
chosen as the focal nuclei for the Si and Fe groups. 
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Corresponding to this reduced set of abundances Q is a reduced set of reac- 
tions, with quasi-equilibrium allowing one to ignore the reactions among the 
members of the QSE groups. Unfortunately, the rates of these remaining reac- 
tions are functions of the full abundance set, Y^, and are not easily expressed 
in terms of the group abundances, Y g . Thus, for each Y g , one must solve for 

— * — * — * — * 

Y K and, by Eq. 27, F , in order to calculate Y Q which is needed to evolve Y g 

— * 

via Eq. 16. Furthermore, Eq. 18 requires the calculation of the Jacobian of Z, 

which can not be calculated directly since Y g cannot be expressed in terms of 
Y g . Instead it is sufficient to use the chain rule, 

QY<3 QfQ QY n , x 

= — — 29 



to calculate the Jacobian. Analytically, the first term of the chain rule product 
is easily calculated from the sums of reaction terms, while the second term 
requires implicit differentiation using Eq. 28. Additional details and compar- 
isons with the full ct-network are demonstrated by Hix et al. [70] (See also 
[71]). 



7 Conclusion 



Computational astrophysics is in the midst of a paradigm shift. For many 
years, one dimensional Lagrangian models, which take advantage of the (near) 
sphericity of self-gravitating fluids, have been the standard for modeling the 
lives, deaths and afterlives of stars. Increasingly, we see these models being 
supplanted by multi-dimensional models which can more accurately include 
the wide variety of phenomena which break spherical symmetry (for example, 
rotation, convection and/or magnetic fields). Present day computational capa- 
bility limits the use of such multi-dimensional models to short-lived episodes 
in the stellar life cycle, but the future will see increasing usage of such models. 
An important limitation of the current multi-dimensional models is the small 
number of nuclear species which are evolved within these models, typically 
less than twenty. This is in contrast to the hundreds or thousands of species 
tracked in the one dimensional models. Removing this limitation is vital if we 
are to continue to compare models to the myriad of observations which re- 
veal the elemental and isotopic composition of the Universe. This will require 
continued improvements in the numerical methods we have discussed. 

In this article, we have reviewed the numerical methods presently used to 
study thermonuclear kinetics in astrophysics and highlighted the directions 
in which we see promise for improvements in speed and accuracy. Nuclear 
compositions are currently evolved within astrophysical simulations via ther- 
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monuclear rate equations or Nuclear Statistical Equilibrium. While NSE so- 
lutions are much more economical, principally because of the much smaller 
number of free parameters that must be evolved, they are applicable to only a 
few situations. The principle difficulty with evolution via rate equation is the 
computational cost, which results primarily from the stiffness of the system 
of nuclear reactions. This extreme stiffness requires implicit solution and has 
thus far generally precluded the use of integration methods that do not rely on 
the Jacobian matrix. Such non-Jacobian methods remain highly sought after 
as a means to reduce the computational cost of nucleosynthesis calculations. 
For Jacobian based integration methods, there remain considerable economies 
to be gained by taking better advantage of the sparse nature of the Jacobian 
and by reuse of the inverted Jacobian. 

Physically motivated approaches can also be extremely useful in reducing the 
computational cost. One such is the use of local equilibria to reduce the size of 
the system of rate equations (and reduce problems with matrix singularity). 
Methods based on local equilibria are applicable to many situations where 
global equilibrium has not been achieved. Though the use of hybrid equilib- 
rium networks is in its infancy, it seems that in many of the situations where 
one would have heretofore used a large network coupled with hydrodynamics, 
the hybrid equilibrium networks provide sufficient accuracy and considerable 
reduction in the computational cost. Ultimately, both physically inspired and 
computationally motivated improvements will be necessary as we strive to 
extend the reach of multi-dimensional simulations. 
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